ged <- read.csv("ged_main_2024.csv") 

library(tidyverse)
library(dplyr)
library(tidyr)
library(lubridate)
library(purrr) 
library(ggplot2)
library(gridExtra) 
library(padr) 
library(plotly)
library(scales)
library(tibble)
library(fftwtools)
library(ecp)

#################
##
## FUNCTIONS
##
#################

subsetting <- function(area, dyad, side, ...) {
  civil <- ged[which(ged$control==area & ged$dyad2==dyad),] %>%
    add_row(date_start = "2016-01-01", date_end = "2016-01-01", .before = 1) %>%
    add_row(date_start = "2019-11-30", date_end = "2019-11-30") %>%
    mutate(across(c(date_start, date_end), ymd)) %>% 
    transmute(id, date = map2(date_start, date_end, seq, by = 'day'), deaths_a, deaths_b, deaths_civilians) %>% 
    mutate(across(c(deaths_a, deaths_b, deaths_civilians), ~ ./lengths(date)))  %>%
    unnest(date) %>%
    group_by(date) %>% 
    summarise(across(c(deaths_a, deaths_b, deaths_civilians), sum, na.rm = TRUE))
  if (side=="a"){
    civil$d <- civil$deaths_a - civil$deaths_civilians
  } else {
    civil$d <- civil$deaths_b - civil$deaths_civilians
  }
  civil <- pad(civil) 
  civil[is.na(civil)] <- 0
  return(civil)
}

divisive.scores <- function(mydata){
  library("ecp") 
  diff <- matrix(c(mydata$d), ncol = 1) 
  ediv <- data.frame(name = paste0("civil.", c("as.op","op.as","sdf.is","is.sdf","as.is","is.as")),
                     val = c(0.25,0.25,0.25,0.75,0.25,1))
  alp = ediv$val[which(ediv$name==deparse(substitute(mydata)))]
  output <- c(e.divisive(diff, R = 499, alpha = alp, min.size=30, sig.lvl = 0.05)$estimates) 
  return(output)
} 

difr <- function(mydata,date,agr,div){ 
  ind <- c() ; dif <- c() ; min.dif <- c() 
  fn <- ecdf(mydata) 
  for (j in 1:length(agr)){ 
    ind <- div[-c(1,length(div))] 
    dif <- as.numeric(agr[j]-date[ind]) 
    if (length(dif[which(dif>=0)])>0) {
      min.dif[j] <- abs(min(dif[which(dif>=0)]))
    } else{
      min.dif[j] <- NA
    }
  } 
  output <- mean(min.dif,na.rm=T)
  return(output)
}

iaaft <- function(x, max.iter=1000){
  N <- length(x)
  try({
    require(fftwtools)
    fft <- function(data, inverse=0) fftw(data=data, inverse=inverse)
  })
  x.sorted <- sort(x)  
  x.amp <- Mod(fft(x)) 
  x.old <- sample(x) 
  for (i in 1:max.iter){
    x.old.ft <- fft(x.old)
    phases <- Arg(x.old.ft)
    x.replace <- fft(x.amp*exp(1i*phases), inverse = TRUE)/N 
    x.final <- x.sorted[rank(x.replace)]
    if (all(x.final==x.old)){
      break
    } else {
      x.old <- x.final
    }
    if (i == max.iter) warning('Reached maximum iterations without convergence!')
  }
  return(Re(x.replace))
}

sig.iaaft <- function(mydata,agr){
  new <- numeric(length(mydata$d)) 
  ediv <- data.frame(name = paste0("civil.", c("as.op","op.as","sdf.is","is.sdf","as.is","is.as")),
                     val = c(0.25,0.25,0.25,0.75,0.25,1))
  alp = ediv$val[which(ediv$name==deparse(substitute(mydata)))]
  min.sur <- numeric(500) 
  for (i in 1:500){
    new <- matrix(c(iaaft(mydata$d)), ncol = 1)
    div.est <- e.divisive(new, R = 499, alpha = alp, min.size=30, sig.lvl = 0.05)$estimates
    min.sur[i] <- difr(new,mydata$date,agr,div.est) 
    print(i)
    timestamp()
  } 
  return(min.sur)
} 


#################
##
## PREPARE DATA
##
#################

## Assad KILLS in Opposition TERRITORY
civil.as.op <- subsetting(area="Opposition", dyad="11973",side="b") 
## Opposition KILLS in Assad TERRITORY
civil.op.as <- subsetting(area="Asad", dyad="11973",side="a") 
####### 

## Assad KILLS IN IS TERRITORY
civil.as.is <- subsetting(area="IS", dyad="14620",side="b")
## IS KILLS IN Assad TERRITORY
civil.is.as <- subsetting(area="Asad", dyad="14620",side="a")
#########

## SDF KILLS IN IS TERRITORY
civil.sdf.is <- subsetting(area="IS", dyad="14637",side="b") 
## IS KILLS IN SDF TERRITORY
civil.is.sdf <- subsetting(area="SDF", dyad="14637",side="a") 
#########

agreements <- c("2016-02-03", "2016-09-10", "2017-01-25", "2017-03-03", "2017-03-31", "2017-07-05", "2017-09-15", "2017-12-22", 
                "2018-01-30", "2018-05-15", "2018-07-31", "2018-11-29", "2019-04-26", "2019-08-02", "2019-11-06") 
agree <- data.frame(dates=agreements,agree=rep(1,length(agreements)))
agree <- agree %>% add_row(dates = "2016-01-01", agree = 0, .before = 1)
agree <- agree %>% add_row(dates = "2019-11-30", agree = 0)
agree$dates <- as.Date(agree$dates, "%Y-%m-%d")
agree <- pad(agree) 
agree[is.na(agree)] <- 0
agreements <- as.Date(agreements, "%Y-%m-%d") 

#################
##
## ANALYSIS
##
#################

set.seed(999)

divisive.as.op <- divisive.scores(civil.as.op)
divisive.op.as <- divisive.scores(civil.op.as)
divisive.as.is <- divisive.scores(civil.as.is)
divisive.is.as <- divisive.scores(civil.is.as)
divisive.sdf.is <- divisive.scores(civil.sdf.is)
divisive.is.sdf <- divisive.scores(civil.is.sdf)

sur.as.op <- sig.iaaft(civil.as.op,agreements) # The next 6 lines takes about two days to run with Mac Pro M3
sur.op.as <- sig.iaaft(civil.op.as,agreements) 
sur.sdf.is <- sig.iaaft(civil.sdf.is,agreements)
sur.is.sdf <- sig.iaaft(civil.is.sdf,agreements)
sur.as.is <- sig.iaaft(civil.as.is,agreements)
sur.is.as <- sig.iaaft(civil.is.as,agreements)

difr_data<-data.frame("asop"= difr(civil.as.op$d,civil.as.op$date,agreements,divisive.as.op),
                      "opas"= difr(civil.op.as$d,civil.op.as$date,agreements,divisive.op.as),
                      "sdfis"=difr(civil.sdf.is$d,civil.sdf.is$date,agreements,divisive.sdf.is),
                      "issdf"=difr(civil.is.sdf$d,civil.is.sdf$date,agreements,divisive.is.sdf),
                      "asis"= difr(civil.as.is$d,civil.as.is$date,agreements,divisive.as.is),
                      "isas"= difr(civil.is.as$d,civil.is.as$date,agreements,divisive.is.as))
